# August 2023
# Statistical analysis of the xylem vulnerability curve and pressure volume curve parameters 
# Statistical analysis of the hydraulic and stomatal safety margin parameters

# Anova and Tukey HSD (posthoc)

# p50
# tlp data
# WP Nov 2020 - very dry
# WP Nov 2022
# 
ov.data$gen_spe <- paste(ov.data$genus,"_",ov.data$species, sep="")
pe.mean <- tapply(ov.data$pe,ov.data$gen_spe,mean,na.rm=T)

# find the water potential difference between P50 for each species and the TLP or WP measurements in the field
# e.g., match the summary.data$tip.label with the tlp.dat$species

# TLP
tlp.dat$ssm <- tlp.dat$TLP.MPa - summary.data[match(tlp.dat$species,summary.data$tip.label),3]
ssm.mean <- tapply(tlp.dat$ssm,tlp.dat$species,mean,na.rm=T)
ssm.se <- tapply(tlp.dat$ssm,tlp.dat$species,se)
# TLP from Pe
tlp.dat$ssm.pe <- tlp.dat$TLP.MPa - pe.mean[match(tlp.dat$species,row.names(pe.mean))]
ssm.pe.mean <- tapply(tlp.dat$ssm.pe,tlp.dat$species,mean,na.rm=T)
ssm.pe.se <- tapply(tlp.dat$ssm.pe,tlp.dat$species,se)

# WP 2020
wp.dry.dat$hsm <- wp.dry.dat$Water_potential_MPa - summary.data[match(wp.dry.dat$species,summary.data$tip.label),3]
hsm.dry.mean <- tapply(wp.dry.dat$hsm,wp.dry.dat$species,mean,na.rm=T)
hsm.dry.se <- tapply(wp.dry.dat$hsm,wp.dry.dat$species,se)
names(wp.dry.dat)
hsm.dry.ss <- tapply(wp.dry.dat$Tree_Ind,wp.dry.dat$species,max)

# WP 2022
wp.dat$gen_spe <- paste(wp.dat$Genus,"_",wp.dat$species,sep="")
wp.dat$hsm <- wp.dat$Water_potential_MPa - summary.data[match(wp.dat$gen_spe,summary.data$tip.label),3]
hsm.mean <- tapply(wp.dat$hsm,wp.dat$gen_spe,mean,na.rm=T)
hsm.se <- tapply(wp.dat$hsm,wp.dat$gen_spe,se)
hsm.ss <- tapply(wp.dat$Individual,wp.dat$gen_spe,max)

# WP 2020
wp.dry.dat$hsm.pe <- wp.dry.dat$Water_potential_MPa - pe.mean[match(wp.dry.dat$species,row.names(pe.mean))]
hsm.pe.dry.mean <- tapply(wp.dry.dat$hsm.pe,wp.dry.dat$species,mean,na.rm=T)
hsm.pe.dry.se <- tapply(wp.dry.dat$hsm.pe,wp.dry.dat$species,se)
# WP 2022
wp.dat$hsm.pe <- wp.dat$Water_potential_MPa - pe.mean[match(wp.dat$gen_spe,row.names(pe.mean))]
hsm.pe.mean <- tapply(wp.dat$hsm.pe,wp.dat$gen_spe,mean,na.rm=T)
hsm.pe.se <- tapply(wp.dat$hsm.pe,wp.dat$gen_spe,se)

# TLP
anova(lm(tlp.dat$TLP.MPa~tlp.dat$species))
x <- aov(lm(tlp.dat$TLP.MPa~tlp.dat$species))
summary(x)
TukeyHSD(x)
# boxplot(tlp.dat$TLP.MPa~tlp.dat$species)
# Result:
# Some variation in TLP between species
# a  Schotia
# ab Polygala
# bc Boscia
# c Pappea,  
# d Searsia, Euclea

# SSM from P50
anova(lm(tlp.dat$ssm~tlp.dat$species))
x <- aov(lm(tlp.dat$ssm~tlp.dat$species))
summary(x)
TukeyHSD(x)
# boxplot(tlp.dat$ssm~tlp.dat$species)
# Results:
# All close stomata prior to experiencing embolism 

# Schotia, Pappea and Polygala have similar, narrow SSMs
# Schotia and Polygala close stomata early, but have least resistant xylem
# Pappea closes stomata at intermediate water potentials, and has relatively vulnerable xylem 

# Boscia, Euclea, Searsia also have similar, large SSMs
# Boscia has intermediate P50, but early stomatal closure 
# Other two species close stomata late, but have very resistant xylem.

# SSM from Pe
anova(lm(tlp.dat$ssm.pe~tlp.dat$species))
x <- aov(lm(tlp.dat$ssm.pe~tlp.dat$species))
summary(x)
TukeyHSD(x)
# boxplot(tlp.dat$ssm.pe~tlp.dat$species)

# HSMs from P50
# During drought
anova(lm(wp.dry.dat$hsm~wp.dry.dat$species))
x <- aov(lm(wp.dry.dat$hsm~wp.dry.dat$species))
summary(x)
TukeyHSD(x)
# par(mfrow=c(2,1))
# boxplot(wp.dry.dat$hsm~wp.dry.dat$species, ylim=c(-2,6))
# abline(h=0,lty=2)
# Results:
# Polygala experienced low water potentials that could have resulted in substantial embolism
# HSM P50 < zero

# Schotia, and Pappea also experienced low water potentials that could have resulted in embolism
# HSM P50 > zero, hsm pe < 0

# Euclea, and Searsia had large HSMs. Unlikely to have experienced stem embolism

# During wetter times
anova(lm(wp.dat$hsm~wp.dat$species))
x <- aov(lm(wp.dat$hsm~wp.dat$species))
summary(x)
TukeyHSD(x)
# boxplot(wp.dat$hsm~wp.dat$species, col="lightblue", ylim=c(-2,6))
# abline(h=0,lty=2)
# Results:
# All species maintained positive hydraulic safety margins in the wetter period
# Polygala and Pappea experienced low water potentials that could have resulted in some embolism
# HSM Pe < zero
# All other species had large HSMs. Unlikely to have experienced substantial stem embolism

# HSMs from Pe
# During drought
anova(lm(wp.dry.dat$hsm.pe~wp.dry.dat$species))
x <- aov(lm(wp.dry.dat$hsm.pe~wp.dry.dat$species))
summary(x)
TukeyHSD(x)
par(mfrow=c(2,1))
# boxplot(wp.dry.dat$hsm.pe~wp.dry.dat$species, ylim=c(-5,5))
# abline(h=0,lty=2)
# Results:
# Polygala experienced low water potentials that could have resulted in substantial embolism
# HSM Pe << zero

# Pappea (and Schotia) also experienced low water potentials that could have resulted in embolism
# hsm pe < 0

# Searsia and Euclea had positive HSMs. Unlikely to have experienced stem embolism

# During wetter times
anova(lm(wp.dat$hsm.pe~wp.dat$species))
x <- aov(lm(wp.dat$hsm.pe~wp.dat$species))
summary(x)
TukeyHSD(x)
# boxplot(wp.dat$hsm.pe~wp.dat$species, col="lightblue", ylim=c(-2,6))
# abline(h=0,lty=2)
# Results:
# All species maintained positive hydraulic safety margins in the wetter period
# Unlikely to have experienced stem embolism